% 1D peridogram analysis for hill profiles - MAIN program 


%   rev 55: Feb 2 2011: added output of detrended profile to csv file
%   rev 54: Jan 26 2011 
%       - removed plotting of the pre avg. pgram 
%       - added 2 doplot args to get_trend for plotting trend functions or residuals
%       - added run configs with undefined variable catching
%       - added line to only plot last of the multi-smoothed pgrams
%   rev 53: Dec 16 2010 - played with using a per-folder parameter file to 
%           read in values to overwrite global settings
%   rev 52: Dec 2  2010 - now %% indicate parameters for tweaking

% reads in elevation profile a tab-separated .txt file (made from ArcGIS's profile), shows
% the (detrended) hill profiles, calculates different types of average peridograms,
% a baseline and 95% (etc.) sigificance lines

% Credits: uses parts of Perron's example.m and David Meko's code

% Run configuration variables: the following variables are supposed to be
% defined in the Run configuration, this will define them if the don't
% exist yet
if exist('fname') == 0 fname = 'prof7c.txt'; end
if exist('doplot_residuals') == 0 doplot_residuals = 1; end
if exist('doplot_trend_funcs') == 0 doplot_trend_funcs = 1; end
if exist('trend_order') == 0 trend_order = 1; end

%% Name of elevation data file 
%fname = uigetfile('*.txt') % open dialog to get text file's name from the user
%fname = 'sample_large_bundle.txt'
%fname = 'prof7_all_even.txt'
%fname = 'prof7c.txt'
%fname = 'profile11.txt'
%fname = 'soi.txt' % used in http://www.colorado.edu/geography/class_homepages/geog_4023_s09/TimeSeriesLab2.html
%fname = 'synth_50m_wave_amp2.0+10m_wave_amp1.0_no_trend.txt' % synthetic profile

foldername = fname(1:length(fname)-4) % folder for output, chop off last 4 letters of filename
%foldername = [foldername '_version2'];  % use this tp add a variation (if you don't want to overwrite existing plots)
disp(sprintf('current output folder is %s', foldername));
mkdir(foldername); % make folder for all plots
fullname = fullfile(foldername, fname(1:end-4)); % full name of files to be written but w/o extension

matrix = textread(fname);
y_cols = matrix(:, 2:2:end); % grab 2., 4., 6., etc. column
x_cols = matrix(:, 1:2:end); 

% % Optional: if t is unevenly spaced, interpolate
% t_even=min(t):dt:max(t);
% y = interp1(t(:),y(:),t_even(:));
% t=t_even(:);
% clear t_even;

num_y_cols = length(y_cols(1,:));
disp(sprintf('read in %d columns of %d samples each', num_y_cols, length(y_cols(:,1) )))

%% Enable plotting of Elevation data
% plot true elevation (prior to detrending), assumes same length of x and y_cols 
%plot_elev_profiles(x_cols, y_cols, foldername); % just for plotting - can be skipped

% get (pre) averages for detrending
y_avg = mean(y_cols, 2); % average elev. profile (2 here means: avg. the 2. dimension, i.e. along columns)
x_avg = mean(x_cols, 2); %
num_samples = length(y_avg); % assumes same length for all profiles!

%% Calculate trend from avg. elevation
%doplot_trend_funcs = 1; % 0 to skip plotting the trend functions
%doplot_residuals = 1; % 0 to skip plotting the residuals
trends = get_trend(x_avg, y_avg, foldername, doplot_trend_funcs, doplot_residuals); % returns an array of trends from the avg. elevation

%% Set order of trend removal
%trend_order =  1 % trend order 
trend = trends(:,trend_order); % copy the column containing the nth order trend
disp(sprintf('detrending with trend order %d', trend_order));

% detrend all columns with choosen trend (with at least 1. order!)
%y_cols_with_trend = y_cols; % save with-trend for later?
for c=1:num_y_cols,
    y_cols(:,c) = y_cols(:,c) - trend;
end
y_avg = y_avg - trend; % detrend avg. elev.
csvwrite([fullname '_detr.csv'],[x_avg y_avg]); % save detrended x/y coords in csv file

%% Sub-setting of elevation (will only use first column, experimental only)
% Mirror-glue it and split into sections
% this should improve unsmoothed pgram as it's avg. a bundle of pgrams of profiles
% i.e. set smoothing_type to 0
if 0 % 
%if num_y_cols == 1 % if we have only 1 column
    l = length(y_cols(:,1));
    %% Number of sub_sets (aprox.)
    num_sub_sets = 100;
    size = floor(l / num_sub_sets);     
    ymg = cat(1,y_cols(:,1), y_cols(end:-1:1,1)); % mirror and glue
    %figure(); plot(ymg); % debug: plot mirrored elevation profile
    c = 1;
    for i=1:size:l,
        y_cut = ymg(i:i+l-1);
        x_cut =[1:length(y_cut)]';
        y_cols(:,c) =  y_cut;
        x_cols(:,c) = x_cut;
        c = c+1;
    end
end

%% Set padding, significance level and win_size for chi-square calc.
dowindow = 1; % tapering - window the data prior to computing FFT, always leave at 1
dopad = 2; % 1 (or more): pad data with zeros to next (next-next) power of 2 
sig_lvl = 0.99; % confidence level used
win_size= 2; % 1 means no smoothing - win_size here is ONLY used for calculating 
            % the degrees of freedom for the chi-square for the
            % theoretical and significnce lines, NOT for the
            % actual smoothing of the pgram, which is done separately later.
                        
if dopad
    % calculate next power of 2 greater than the number of samples 
    % actual padding is done inside fft1d()
    extra = dopad-1;
    num_bins = 2.^((ceil(log(num_samples)/log(2)))+ extra);
else % no  padding, 
    num_bins = num_samples;
    num_bins = floor(num_bins / 2) * 2; % make it an even number    
end


% Start new plot
figure 

% calc. peridograms and signif lines of individual elevations
for c=1:length(y_cols(1,:)),
    y = y_cols(:, c); % current column
    %t=1:length(y); % "time" (dt=1) % same as x?
    t = x_avg;   % instead of a sample step size of 1, use real distances
    [P f] = fft1d_pgram(t,y,dopad,dowindow);% get pgram and f for this profile, f will be identical for all
    [fft_theor signif lag1]=fft1d_theor_signif(P,f,t,y,sig_lvl,win_size);
    P_indiv(:, c) = P; % store individual results as column c in a 2D array
    fft_theor_indiv(:, c) = fft_theor;
    fft_signif(:, c) = signif;
end

%% Make and plot "pre' peridogram gram from the averaged elevation
t=x_avg;
[P_avg_pre f] = fft1d_pgram(t,y_avg,dopad,dowindow);% pgram of avg. elevations (y_avg)
sig_lvl = 0.99;
[fft_theor_avg_pre signif_avg_pre lag1]=fft1d_theor_signif(P_avg_pre,f,t,y,sig_lvl,win_size);
% calculate the fraction of the spectrum that exceeds the significance
% level. This should be approximately equal to 1-siglvl, or 0.05 if siglvl=0.95
sig_frac_pre = sum(P_avg_pre>signif_avg_pre)/length(P_avg_pre); % fraction exceeding significance

%loglog(f,fft_theor_avg_pre,'-m') % seems to always == fft_theor_avg_post, so I'm omitting it
%loglog(f,signif_avg_pre,'-.m')

%%  Plot "post" peridogram (average of indiv. pgrams)
t=x_avg; 
P_avg_post = mean(P_indiv, 2); 
fft_theor_avg_post = mean(fft_theor, 2);
signif_avg_post = mean(fft_signif, 2);
sig_frac_post = sum(P_avg_pre>signif_avg_post)/length(P_avg_post); % fraction exceeding significance

h_P_avg_post = loglog(f,P_avg_post,'xk','LineWidth',1, 'color', [1 1 1]); % keep plot handle for later
hold on

theor_color = [0.7 0.3 0.6]; % color for theoretical and baseline
h_fft_theor_avg_post = loglog(f,fft_theor_avg_post,'-', 'color', theor_color);
h_signif_avg_post = loglog(f,signif_avg_post,'--', 'color', theor_color);
%h_P_avg_pre = loglog(f,P_avg_pre,'-','LineWidth',0.5, 'color', [0.5 0.5 0.5]);  % pre avg pgram 

%% Assemble a legend from plots
legend_strings = {'peridogram (post)',...  % strings for each plot, must use {} here!
    sprintf('theor. spectr. (win.sz=%d)', win_size), ... 
    sprintf('%0.2f%% signif. (%0.3f abv.sigf.)', sig_lvl, sig_frac_post), ...
    %'peridogram (pre)'
    };

% combine plot handles into a "package"
%plot_handles = [h_P_avg_post, h_fft_theor_avg_post, h_signif_avg_post, h_P_avg_pre]; 
plot_handles = [h_P_avg_post, h_fft_theor_avg_post, h_signif_avg_post]; % removed h_P_avg_pre 

legend(plot_handles, legend_strings) % make a legend from the handles and their stings
    %'Location','SouthWest') %  bottom left
    %'Location','Best') % set location to adjust to empty space
    
    
%% Set smoothing type (none, single or multi-pass) 
filter_name='none'; wins = [0]; % default: no smoothing
hold all  % difference to on is that all cycles through the colors

% note:  wins contains the sizes of the smoothing window (filter width), MUST be odd numbers!
%        the window type used (rectangle, daniell), is set inside apply_filters() or re_apply_filters()    
smoothing_type = 2; % 0=no smoothing, 1=single-pass 2=multi-pass
%wins = [3,5,7,11,15, 19, 21]; % vector (list) of filter window widths (smoothing type = 1)
wins = [3,3,3,3,3,3,3,3,3,3]; % vector (list) of (succcessivly) filtered window widths (smoothing type = 2)
%wins = [21]; % only a single smoothing pass

if smoothing_type == 0 % No smoothing
    for c=1:length(y_cols(1,:)),
        %loglog(f,P_indiv(:, c), 'xk','MarkerSize',3, 'LineWidth',0.001) % debug: plot all indiv. pgrams
        %loglog(f,P_indiv(:, c), 'LineWidth',0.001)
    end
    loglog(f,P_avg_post,'k','LineWidth',1.5); % plot post_avg pgram
end

% set the colormap for smoothed pgrams
if smoothing_type >= 1
    colormap(hsv(length(wins))) % use the (hsv, summer, winter, jet, etc. colormap for smoothed pgrams
    cm = colormap;  % type colormapeditor to play with colormaps
    set(gca, 'ColorOrder', cm); % makes
    hold all
    colorbar('location', 'SouthOutside',...
        'Units','Normalized',...
        'Position', [0, 0.0 1 0.05],... 
        'OuterPosition',[0, -0.0001 1 0.03],...
        'Xtick', [])
end

%%  Set filter type for single pass smoothing inside apply()
% each element of wins is used to smooth the raw pgram and plot the result
if smoothing_type == 1
    [P_sm filter_name] = apply_filters(P_avg_post, wins); % 1 or many smoothed pgrams (array)
    loglog(f, P_sm(:,end), 'LineWidth',0.5); % plot the last of the single-pass smoothed pgrams
    disp(wins);
end

%% Set filter type for multi-pass smoothing inside re_apply()
% here the raw pgram is smoothed and the re-used for the 2. (3., 4., ...) pass
if smoothing_type == 2 
    [P_sm_r filter_name] = re_apply_filters(P_avg_post, wins);
    %loglog(f, P_sm_r, 'LineWidth',0.5); % plot pgrams as they are succcessivly smoothed
    loglog(f, P_sm_r(:,end),'r-', 'LineWidth',0.5); % plot only last smoothed pgrams (in red)
    disp(wins);
end

%% Set plot parameters 
set(gcf,'Position',[0 0 1200 800])% window size: upper left corner position (x,y), width,  height

% limit view?
xlim([1/1000 1/40])  % must set xlim in fractions of dist here - not at the  wavelength axis later!
%xlim([1/100 1/5])
%xlim([1/500 1/50]) 
%ylim([min(P_avg_post) max(P_avg_post)]) % sloppy way of setting of y range

axf = gca;  % get x axis for frequency
frange = get(axf,'xlim');
wrange = 1./fliplr(frange);
set(axf,'box','off','tickdir','in')

% add another x-axis axw on top for the wavelength
axw=axes('Position',get(axf,'Position'),...
         'XAxisLocation','top',...
         'YAxisLocation','right',...
         'Color','none',...
         'XColor','k','YColor','k',...
         'xlim',wrange,'xdir','reverse',...
         'ytick',[],'xscale','log');
set(axw,'box','off','tickdir','in')

%set(axw,'XTick', [0:10:150]) % custom, even, tick marks - adjust as needed
ticks = horzcat([1:1:15],[20:10:90],[100:25:200],[250],[300:100:length(f)]); % or use a list of ticks
%ticks = horzcat([1:1:19], [20:2:29], [30:5:50]); % or use a list of ticks
%ticks = horzcat([1:1:19], [20:2:29], [30:5:49], [50:10:100]); %
set(axw,'XTick', ticks)
set(axw, 'XMinorTick', 'on');

% these axis label positions are in normalized (0-1) units
set(get(axf,'xlabel'),'string','Frequency (m^{-1})', 'Units', 'Normalized', 'Position', [0.5 -0.02]) 
set(get(axw,'xlabel'),'string','Wavelength (m)', 'Units', 'Normalized', 'Position', [0.5 1.03])
set(get(axf,'ylabel'),'string','Mean squared amplitude (m^2)')
set(axw, 'XGrid', 'on'); % turn on x axis grid for wavelength

% textbox with parameters, text must be in {} without commas!
text(0.99,0.01, ... 
    {
     sprintf('Filter=%s, type=%d, passes:', filter_name, smoothing_type)...
     sprintf('%d ',  wins)... % list of fiter window widths
     sprintf('order of removed trend =%d',trend_order)...
     sprintf('dopad?=%d,   %d bins,   %d samples', dopad, num_bins, num_samples)...
     %sprintf('dowindow=%d', dowindow)...
     %sprintf('win_size=%d', win_size)...
    },...  
    'Units', 'Normalized', 'VerticalAlignment', 'bottom','HorizontalAlignment', 'right')

title_string = ['Data: ' fname];
title(title_string, 'Units', 'Normalized', ...
    'Fontsize',14, 'Interpreter','none');



if 1 % print into files 
    %set(gcf,'color','none'); % set background color transparent (works on mac for pdfs, may not work on PC) 
    
    %set(gcf, 'PaperPositionMode', 'manual');
    set(gcf, 'PaperUnits', 'inches');
    set(gcf, 'PaperType', 'usletter');
    set(gcf, 'PaperPosition', [0 0 11.5 8]);
             
    % for file formats: http://www.mathworks.com/help/techdoc/ref/saveas.html
    set(gcf, 'PaperOrientation', 'landscape'); % pdf is chopped off otherwise 
    saveas(gcf, [fullname '.pdf']);
    set(gcf, 'PaperOrientation', 'portrait'); 
    %saveas(gcf, [fullname '.emf']); % won't work on Mac
    saveas(gcf, [fullname '.eps']); % encapsulatd postscript (vector)
    saveas(gcf, [fullname '.png']);
    
end
hold off



% clear  % clear workspace

